
use $datapath/clean/main_data.dta, clear 


************************************************
* merge with income info 

merge m:1 fipsplac fipsst using "$datapath/intermediate/income_place_data.dta"
	drop if _merge==2
	drop _merge 

	
************************************************
* Create variable labels: 

	label var Qp "phys"
	label var Np "phys holdings"
	label var Qe "ebook"
	label var Ne "ebook holdings"


** three income levels 
	xtile ino = medinc, nq(3)


*****************************		

**********************************************************************

**********************************************************************
** Reshape data to run one regression for all

	keep Qp Qe Np Ne year lno ino
	gen num=_n 
	reshape long Q, i(Np Ne year lno num) j(type) s 

	gen de = type=="e"
	egen lfno = group(lno de)
	egen yfno = group(year de)




*******************************************
* Income terciles (higher = more)

	egen lfino = group(lno de ino)
	egen yfino = group(year de ino)

areg Q i.de#i.ino#c.Np i.de#i.ino#c.Ne i.yfno, absorb(lfno)  vce(cluster lno)

** thetas: 	
	nlcom ((2/3)*_b[0.de#1.ino#c.Np] - _b[0.de#1.ino#c.Ne]) / (_b[1.de#1.ino#c.Ne] - (2/3)*_b[1.de#1.ino#c.Np])
	nlcom ((2/3)*_b[0.de#2.ino#c.Np] - _b[0.de#2.ino#c.Ne]) / (_b[1.de#2.ino#c.Ne] - (2/3)*_b[1.de#2.ino#c.Np])
	nlcom ((2/3)*_b[0.de#3.ino#c.Np] - _b[0.de#3.ino#c.Ne]) / (_b[1.de#3.ino#c.Ne] - (2/3)*_b[1.de#3.ino#c.Np])
	
